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Abstract. We study the dynamics of density fluctuations in the steady state of a 
non-equilibrium system, the Zero-Range Process on a ring lattice. Measuring the 
time series of the total number of particles in a segment of the lattice, we find 
remarkable structures in the associated power spectra, namely, two distinct components 
of damped-oscillations. The essential origin of both components is shown in a simple 
pedagogical model. Using a more sophisticated theory, with an effective drift-diffusion 
equation governing the stochastic evolution of the local particle density, we provide 
reasonably good fits to the simulation results. The effects of altering various parameters 
are explored in detail. Avenues for improving this theory and deeper understanding of 
the role of particle interactions are indicated. 
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1. Introduction 

For systems in tliermal equilibrium, Boltzmann and Gibbs provided a sound framework 
which forms part of text-book material nowadays. Specifically, the probability 
distribution (for finding the system in any configuration) is given simply by the 
Boltzmann factor. In contrast, little is known in general about systems driven into 
non-equilibrium steady states, specified via, say, a set of transition rates that violate 
detailed balance. To explore this vast unknown, it is natural to focus on simple solvable 
models, with the hope of gaining insight for formulating a general framework for non- 
equilibrium statistical mechanics. One paradigmatic model is the Zero-Range Process 
(ZRP) [11 [2], in which particles hop from one site to the next on a one- dimensional 
periodic lattice (i.e., a ring of L sites), with a rate that depends solely on n, the number 
of particles in the originating site. 

The ZRP distinguishes itself in at least two ways. Its steady-state probability 
distribution not only takes a factorised form, which can be expressed succinctly in terms 
of the hopping rate u{n) (often easing the computation of various averages). It also 
exhibits condensation transitions, even in low dimensions which would not be expected 
in an equilibrium system without long-range interactions. In addition to its utility in 
the fundamental study of non-equilibrium processes, where it has been employed to 
develop a general criterion for phase separation in one- dimensional driven systems [3] 
among other things, the ZRP has also found success as a minimal model for various 
real systems including vehicular traffic [1], compartmentalised granular gases [5] and gel 
electrophoresis [6]. 

Though much is known about the ZRP, there is a simple and natural, but so far 
unknown, question we may ask. While the total number of particles on the (finite, 
periodic) lattice is fixed as the system evolves, the number in a subsection of, say, £ 
sites - denoted by A^^ (t) here - is a quantity that fluctuates in time t. In the steady 
state, its time average {Ni(t)) is of course a constant. Nevertheless, its average power 



autocorrelation {Ne (t) Ni {t')) . A recent study [7] reported the presence of interesting 
oscillations in the power spectra for another simple paradigmatic non-equilibrium model, 
the open Totally Asymmetric Simple Exclusion Process (TASEP) [HIHIE]- In this work, 
we carry out an extensive investigation of the ZRP, with a variety of rates and a range 
of ring and subsection sizes (L and i). In addition to the oscillations discovered in 
TASEP on an open lattice [7j, we found two, "complementary" sets of oscillations, one 
controlled by L and the other, by i. If we keep the subsection length flxed and let the 
ring size go to infinity, we would have the equivalent of an open lattice (with appropriate 
input and outfiux rates to ensure a finite and non- vanishing {Ni (t))). In this sense, our 
considerations for a closed, periodic lattice can be regarded as inclusive of open lattices 
as well. 

The dynamics of density fiuctuations in the ZRP has been investigated recently by 
Gupta, et. al. [lO], who focused on the variance of the integrated current through a single 




non-trivial and provides information on the 
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site as a function of time. This quantity also displays damped oscillations, but in the time 
domain. The period is proportional to L/v, with v being the velocity associated with 
a fluctuation. Being in the frequency domain, our study complements the earlier work. 
One of the components of our I {oj) oscillates with period v / and undoubtedly can be 
traced to the same origin. However, in [10] , the damping of the oscillations is not the 
main focus and so, diffusive/ dispersive properties in the system were not considered. By 
contrast, we will analyse such behaviour and find that they pose the most challenge for 
theoretical understanding. Further, by considering a quantity associated with multiple 
sites (£ > 1), we hope to extract more information about the dynamics of fluctuations. 

The presence oi v/L may lead some to dismiss these phenomena as "nothing but 
finite size effects." However, to understand these effects is very important, since many 
such models of non-equilibrium transport are believed to applicable to physical systems 
of relatively small Ls. Thus, in the examples given above, L < 10, 000 for vehicular 
traffic and gel electrophoresis. The other popular model of non-equilibrium transport, 
TASEP, was first introduced as a possible model for protein synthesis [H], where L 
rarely exceeds 1000. By contrast, in traditional macroscopic systems (e.g., in solid 
state physics), we generally have in mind ~ (10^)^. Consequently, the phenomena 
discussed here are not merely of theoretical interest, but should be physically observable. 

The paper is organised as follows. Details of the model and its dynamics will be the 
focus of the next section. Results of simulations, theoretical analysis and comparisons 
are the themes of the following three sections. We end with a summary and outlook for 
future research. 

2. Model specification and simulation details 

The ZRP studied here comprises a one-dimensional lattice of L sites with periodic 
boundary conditions. A total number of particles is placed on the lattice, with 
no restrictions on the occupation of each site, so that a configuration of the system 
is completely specified by the set {n(x)}, i.e., n (x) particles being on site x (with 
X = 0,...,L — 1). This system evolves through particles hopping from site to site. 
Particles hop to the rightmost adjacent site with a rate that depends only on n: u{n). 
That this hopping rate is independent of the occupation in any other site in the system 
gives rise to the ZR part of ZRP. A basic diagram of the system is shown in Figure [H 
Ours is perhaps the simplest ZRP, with the particle landing only in the nearest 
neighbouring site. Many more complex cases exist, e.g., having more than one particle 
move, landing them in various sites, and inhomogeneous hopping rates: u{n,x). For a 
recent review see [2j. 

Obviously, n{x) = N is a. constant in time. The most naive expectation of our 
simple ZRP is that, after the system settled into a (non-equilibrium) stationary state, 
the average occupation is homogeneous: {n{x)) = p = N/L. Part of the general interest 
in the ZRP arises in the existence of a phase transition, when p exceeds a critical value pc, 
from such a homogeneous state to an inhomogeneous "condensed" state. In the latter. 
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Figure 1. General diagram of the system. Particles on a lattice hop to the rightmost 
adjacent site with a rate that depends on the number of particles at the departure 
site, u{n). The lattice has L sites and periodic boundary conditions. The quantity of 
interest is the total number of particles in a segment of £ sites of the lattice. 



all but one site is occupied by an average of pc particles, with the excess (A^ — pcL) in 
just a single site. Reminiscent of Bose-Einstein condensation, this macroscopic fraction 
is referred to as the condensate. Unlike the Bose-Einstein case, translational invariance 
is spontaneously broken and the condensate can reside in any site. Indeed, in a finite 
system, it does disappear from one site and reappear in another on some interesting 
time scale jT2j. 

Of the infinitely many functions we can choose for u (n) , it is the n ^ oo asymptotic 
properties that control the existence of a transition. A typical rate is u{n) = 1 + b/n, 
which allows a condensation transition provided b > 2 [13]. In this paper, we consider 
several rate functions, including this one, a constant hop rate, and u{n) oc n (which 
corresponds to having noninter acting particles in the lattice). Despite the interest in 
condensation and phase transitions, here we will focus on the homogeneous phase of the 
system, in which many remarkable features already appear. In future studies, we plan 
to further investigate the power spectra of systems with a condensate, as well as lattices 
with more interesting topologies [H]. 

The system is studied using a simple Monte Carlo algorithm. In each Monte-Carlo 
Step (MCS), we make L attempts to move a particle. In an attempt, a site is selected at 
random. Provided it contains n (> 0) particles, one is moved to the rightmost adjacent 
site with probability 'yu{n), where 7 is a normalisation factor, < 7 < 1/ max 'u(n). Of 
course, a particle leaving site L — 1 is moved to site 0. Simulations were typically run 
for 8 X 10^ MCS. 

Starting from random initial conditions, we discard the first 10^ MCS for the system 
to come to the steady state. Subsequently, we focus on a fixed segment of length £ sites 
and record its total occupation, 

i-i 

N,{t)=J2^{x,t) , (1) 

x=0 

every 10 MCS. Thus, each of our time series consists of 7 x 10^ data points, which we 
regard as 53 samples of 131072 points, i.e., t = 0, . . . , T — 1 with T = 2^"^. We choose a 
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power of 2 so that fast Fourier transform routines can be exploited to compute 

N^{uj) = ^e-''^^Ni{t) ; uj = 2'nm/T. (2) 

t 

Finally, we average over the 53 samples to arrive at the power spectrum: 




A typical system size of L = 10000 sites was used with varying segment sizes, but 
most frequently ^ = 1000. Mostly, we use the Mersenne Twister [15] random number 
generator. To rule out systematic errors from these sources, we have also used other 
generators (e.g., drand48, ran2, /dev/urandom), as well as various data segment sizings 
and output intervals. 

3. Simulation results 

A typical trace of Ni{t) is presented in Figure [21 showing both an entire sample and 
a magnified view of a section of this sample. It appears that there is some oscillatory 
behaviour, but it is difficult to distinguish from random fluctuations; taking averaged 
power-spectra measurements reveals much more in terms of structure. A typical power- 
spectra measurement is shown in Figure [3] for a system with L = 32000, i = 1000 and 
hop rate u{n) = 1 + 4/n. It is a log-log plot and the power-spectrum is plotted against 
the index, m, which is related to the frequency, u, through the relation u = 27rm/T. The 
averaged spectrum displays several prominent features. Two distinct damped-oscillation 
components can clearly be seen: one at low m and the other at higher m. The former 
consists of a series of sharp peaks, with the first peak at m = 14. Note the positive 
curvature at low m, a feature notably absent from previously observed power spectra in 
open systems [7]. The higher-m oscillations are more subtle, being obscured partly by 
the other component. Their character is different, resembling those observed in [7], i.e., 
"dips" over a smooth background. The first "dip" can be seen at m ^ 330, where the 
low-m oscillations are effectively damped out. For large m the spectrum tends to m~^, 
characteristic of white noise that might be expected for frequencies associated with a 
microscopic time-scale. 

The effects of various parameter variations on the structure shown in the power 
spectra were studied numerically. These results are discussed in some detail below. 

Density — Examining the changes elicited in the power spectrum when varying 
the density in the system it becomes clear that increasing the overall density lessens and 
eventually the low-m damped oscillating component. Figure HI However, the same effect 
is not conclusively observed for the higher-m damped oscillations. In fact, simulations 
at higher values of b indicate that this structure remains well into the condensed region, 
but at high enough density the oscillations will no longer be apparent. It is also notable 
that the removal of the low-m oscillations coincides with the onset of condensation in the 
system. For the finite systems studied here condensation does not occur with a sharp 
transition at the theoretical critical density (in this case pc = 0.5) but rather a region 



Power Spectra in a ZRP 



6 



(a) 

350 




150' ' ' ' ' ' ' ' 

2 4 6 8 10 12 

Time(IOMCS) ^^q4 

(b) 



350 




150' ' ' ' ' ' 

2000 4000 6000 8000 10000 

Time (10 MCS) 

Figure 2. A typical trace of the measurement of Ni{t), the total number of particles 
in a section of 1000 sites in a lattice of 10000, from a simulation of a system with 6 = 4 
at a density of p = 0.25. (a) One sample set of 131072 data points, (b) A portion of 
10000 of this set. 
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Figure 3. A typical result of a power-spectra measurement. Taken from a system 
of L = 32000 sites with a segment of £ = 1000, a density of p = 0.25 and hop rate 
u{n) ~ 1 + 4/n. The two damped-oscillation components in the power spectrum can 
clearly be discerned. The first peak of the low-m component is at m = 14 and while 
the location of the first peak of the higher-m component is obscured by the oscillation 
of the other, the second peak is at m « 650. 
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Figure 4. Results displaying the effect of varying the density with the standard hop 
rate. Data shown are for a system of L = 10000 sites with a segment of ^ = 1000 and 
hop rate 1 + 4/n. The density varies from p = 0.25 to p = 0.55 in increments of 0.1. 
Note that the data have been scaled on the y-axis for comparison and so the units of 
this axis are arbitrary and it is a log- log plot. 

of unstable wandering condensates which become more and more stable as the density 
is increased. In the data shown in Figure H] there is a moderately-stable wandering 
condensate for the overall density p = 0.55. 

Segment Size — In Figure [H results are shown for variation in the size of the 
segment {£) in the system, keeping the total number of sites (L) and particles (N) 
fixed. Here it can be seen that increasing the size of the segment causes the higher-m 
oscillation component to move to lower and lower m where its interference with the other 
component is increasingly apparent. Thus, it is clear that the structure of the spectrum 
at high m is due to the size of the segment. Note that for a segment of 5000 sites in the 
ring of 10000, the higher-m component interferes constructively and destructively with 
alternating peaks of the low-m component. 

Lattice Size — The effect of changing L, the lattice size, was also investigated 
and results are shown in Figure [61 It is clear that at fixed segment size and particle 
density, increasing the lattice size changes the low-m damped oscillating component 
but has little effect on the higher-m component. In conjunction with the results for 
segment size, this leads to the conclusion that the higher-m component is controlled by 
the segment size and that the low-m component is controlled by the size of the lattice. 

Parameter b — Remaining with the standard hop rate, the effect of varying the 
parameter b (which can be thought of as controlling the strength of the interaction) 
was investigated. Results for this are shown in Figure [71 It is clear that changing this 
has an effect on the location of the peaks in both the low-m and higher-m components. 
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Figure 5. Results displaying the effect of varying the segment size with the standard 
hop rate. Data shown are for a system with L = 10000 sites, a density of p — 0.25 and 
hop rate 1 + 4/n. The segment size varies from 1000 to 5000 in increments of 1000. 
Note that the data have been scaled for the purposes of comparison, so the units of 
the y-axis are arbitrary and it is a log- log plot. 




Figure 6. Results displaying the effect of varying the lattice size with the standard 
hop rate. Data shown are for a system with a segment of £ = 1000 sites, a density of 
p = 0.25 and the hop rate 1 + 4/n. The lattice sizes shown are 6000, 7000 and 8000. 
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Figure 7. Results displaying the effect of varying the parameter h in the standard hop 
rate u{n) = 1 + h/n. Data shown are for a system with a lattice of L = 10000 sites, a 
segment size oi I — 1000 and density p = 0.25, with h values 2, 3 and 4. The units of 
the y-axis are arbitrary as the data have been scaled so that they may be compared 
on the same graph and it is a log-log plot. 

Also, although not immediately apparent from the data presented, it has an effect on 
the height of the peaks and the number of clearly resolved peaks. 

Hop Rate Form — Moving on from the standard hop rate, a comparison of this 
with constant and non-interacting hop rates is shown in Figure [HI From these results it 
appears that the two damped oscillation components are most clearly seen in the case 
of the non-interacting hop rate and least clearly seen in the standard hop rate. 

The effect of varying the segment and lattice sizes has much the same effect on 
the constant hop rate and noninteracting cases as it did with the standard hop rate. 
However, changing the density continues to affect the power spectra in the case of the 
constant hop rate, but not for the noninteracting case. In the latter, changing the density 
merely changes magnitude of the power-spectrum, as shown in Figure [H This suggests 
that damped oscillations are universal phenomena in particle-transport systems, though 
inter-particle interactions will affect the detailed properties. Indeed, we will show in the 
next section that a drift-diffusion type interpretation is quite successful in describing 
this phenomenon and that the interaction of the particles controls the values of the drift 
and the diffusion. 

4. Theoretical understanding 



Before we present a theory based on the Langevin equation, let us consider a simple toy 
model, from which we can gain some insight into the origins of both types of oscillations. 
In this toy model, a single particle moves on a ring of length L with uniform velocity 
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u(n) = 1 + 4/n 

u(n) = 1 
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Figure 8. Results comparing power spectra from three different hop rates, from top 
to bottom: standard u{n) — 1 + 4/n; constant = 1; and non-interacting u{n) oc n. 
The units on the y-axis are arbitrary as the data have been scaled for easy comparison 
and it is a log-log plot. The data were taken from a system with L = 1000 sites, a 
segment size of ^ = 100, a density of p = 0.1 and the simulations were all run with the 
same normalisation, 7. 
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Figure 9. Results displaying the effect of varying the density in the system with the 
noninteracting hop rate: u{n) oc n. Data shown are for a system of size L = 500 sites 
with a segment of £ = 50 and densities 0.1, 1.0 and 4.0. 
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Vq, and we seek the power spectrum of (t), the number of particles in a segment of 
length i. (For simplicity, assume continuous space-time.) Of course, {t) is trivially a 
series of step functions, such that 

dtNe (t) ^Y,[S{t- f,TL) -5{t-f,TL- n)] (4) 

where 

tl = L/vq and = ^/vq (5) 

is the time it takes to traverse the ring and the segment, respectively. Taken over a 
period of M traversals {t e [0, Mtl), to be precise), {uj) = J e'"^^N(, (t) is given by 

iuN, {u) = J2 [1 - e-'^^^] = ^ _ [1 - e--^] . (6) 



Thus, 



i"toy {(^) = 



sin {ujtlM/2) 
sin {utl/I) 



"' ^ '"sin {ujTi/2y ^ 



(7) 



For M sufficiently large, the first of these factors gives a scries of "spikes" when uo is 
an integer multiple of lnv^jL. Thus, these peaks arc controlled by the ring size L. 
By contrast the second factor displays an oscillation over a smooth background (cu^^), 
noticeably with zeros at integer multiples of 2t^vq/1. We see that these are governed by 
the segment size I. As a toy model, it also serves as a pedagogical tool. What we have 
here is the temporal version of the diffraction pattern from a large (M) array of slits of 
width I, spaced a distance L apart. Of course, once we add dispersion (velocity here; 
wavelength in diffraction) , we will see both smoothing of the peaks and in- filling of the 
zeros. 

With this insight, we turn to the power spectra here, which can be reasonably 
well understood through a Langevin equation for the local particle density (continuous 
variable) on discrete space-time: p{x,t). To connect with the above section, we may 
think of p{x,t) as a kind of coarse-grain average of n{x,t). The starting point is a 
discrete continuity equation 

p (x, t + 1) - p (x, t) = J{x - 1, t) - J(x, t) , (8) 

with J{x^ t) being the net local current from site a; to a; + 1 at time t. Clearly, it is 
controlled by the hop rate u{x,t). Now, as we are considering power spectra for a; > 0, 
we need only account for the deviations of this density from the mean, i.e., 

<f = p{x,t)-p;x = 0,l,...,L-l;t = 0,l,...,T-l. (9) 

Except for the condensed phase, p is just the global density of particles, N/L; otherwise, 
it is pc (apart from the site with the condensate). The strategy is, for systems far from 
criticality, these deviations should be small and their essentials can be understood via 
an approximate Langevin equation that is linear in ip {x,t). 
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Following the standard route, we recognise the deterministic part of J as a function 
of p and expand that to first order in : 

Jdet{x, t) = Jdct(p + <^{x, t)) = Jdct(p) + Jdet(p)v5(a;, t) + ... . (10) 

Defining 

V = 4et(p) (11) 

and adding the noisy part of the current, rj (x, t), we have 

(p{x, t + 1) - (p{x, t) = V [(p{x -l,t) - (p{x, t)] + ■r]{x - l,t) - ri{x, t) . (12) 
The noise is assumed to be uncorrelated Gaussians, so that 

{r]) = 0; {riix, t)r]ix', t')) = AS,,,,St,t' • (13) 

Here, A is a measure of the strength of the noise, which we regard as a phenomenological 
parameter. Before proceeding to the solution, let us emphasise that this somewhat 
unusual form of the drift- diffusion equation is a signature of the ZRP. In many other 
driven diffusive systems [16], the current from x to a; + 1 would depend on both p {x, t) 
and p (x + 1, t), so that an additional term involving (f{x+l,t) will appear on the right 
hand side of equation (fT2l) . By contrast, in the ZRP, the jump rates (from x to x + 1) 
depend only on the occupation at x. 

In this linear approximation, our Langevin equation (|T2l) can be easily solved by 
Fourier methods. Defining 

x,t 

where k = 27rj/L, uj = 2iTm/T (j = 0, 1, 2, . . . , L — 1 and m = 0, 1, 2, . . . , T — 1), we 
find the solution easily: 

Vik, ^) = — r—k T^vik, ^) . (15) 

Note that, if we keep terms to lowest (relevant) order in k and u, the propagator assumes 
the familiar form, 

—ik 



iuj + ivk + Dk^ ' 

of a drift-diffusion equation with conserved noise: dtp = DV^p — vV p — Vrj. For us, the 
zero-range aspect of the ZRP imposes a relation between v and the diffusion "constant" 
D. 

Now, our focus here is the number of particles in a segment and so we consider 

i-i 

N,{t)=pi + J2vM . (17) 

x=0 

Carrying out the sum over x of e'^''^, we find (for u > 0) 



NAuj) = > — : — -rfjik, Uj) 
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from which we can compute the power spectrum via equation ([3]). Before comparing such 
a resuh to data, we recall that, in the previous study [7], the diffusion coefficient seems 
to be seriously renormalised by interactions. Anticipating the same behaviour here, we 
relabel the coefficient of our [1 — cos(A;)] term (i.e., terms even in k) in the propagator as 
— the effective diffusion constant — and regard it as a phenomenological parameter. 
The resulting power spectrum, after performing the average over the noise (equation [13]), 
is 



LT ^ {cos(cj) - 1 - Deff(cos(A;) - 1)}^ + {sin(cj) + v ^m(k)f 

As we will discuss below, to fit the data well, we must choose values of both A and Dgfr 
to be far from those naively derived above. By contrast, we will see that the "bare" 
value of V (i.e., J'detip)^ with Jdet(p) being the known current density relationship of the 
ZRP [131 12]) fits the data quite well. At this stage, there is no good explanation for 
why A and are significantly "renormalised," while v seems to be "unscathed." Our 
conjecture is that Galilean invariance imposes a Ward identity, as in the case of the 
driven lattice gas [IT]. This is an avenue which we plan to pursue in the future. 

Returning our attention to equation (fT9l) . we see that it predicts the locations of 
the first set of peaks (low u ones resulting from the sojourn time of a fluctuation around 
the entire lattice, L) of the power spectra measured from simulation. For example, if we 
consider small (positive) uj, we find peaks in I {uj) whenever sin(u;) +vsm{k) vanishes, 
i.e., u ~ 2TTjv/L {j = 1,2,3, .. . associated with the k = L — j terms). In a similar 
way, cos{ki) introduces oscillations on the scale controlled by the segment length, i. In 
particular, the factor 1 — cos(H) suppresses the j^^ peak if £ is a unit fraction of L, i.e., 
i = L/j. This behaviour, which is reminiscent of interference, is most pronounced in 
the case of j = 2 (top curve L = 10000; £ = 5000) in Figure [5], where the second peak 
in the other curves is clearly "missing." A less prominent "interference" can be seen 
for the third peak in the middle curve (L = 10000; £ = 3000; j ^3). If £ < L, then 
such effects will be noticeable only at cj's that are large compared to those among the 
peaks. In this regime, the damping is so severe that the peaks dissolve into a smooth 
background and the effects of this factor appear as "dips." 

Turning to the specifics of the ZRP, the current (in the thermodynamic limit) is 
equal to the fugacity z [131 12], so that v will be given by dz/dp. In the case of the 
standard hop rate, u{n) = 1 + b/n, the p-z relation is 

z 2Fi{2,2-2 + b;z) 
^ l + b2Fl{l,l;l + b;z) ■ ^ ' 

so that the velocity is 

y = (1 + 6)2^1(1, 1; 1 + 6; z)/ 



17(0 oo^h 4^ I7f',',',^h ^ ^ 2FH2,2;2 + b;z) 

2Fi(2, 2;2 + b;z) + 7^—^2^1(3, 3; 3 + 6; z) - 



2 + b' ' ' ' ^ 1 + 62^1(1,1; 1 + 6; ^) 



(21) 
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Similarly, it is easy to show that, in the other cases, 

1 for M oc n, i.e., non- interacting particles , , 

(1 + p)^^ for u = const. 

With these analytic results, we are in a position to examine how well this theory fits 
the simulation data. 



5. Comparison between analytic and simulation results 

There is generally good agreement between the simple theory presented above and 
results from simulation, i.e., appropriate choices of the parameters A and D^s produce 
reasonable fits. Not surprisingly, the agreement is best for the noninteracting case. In 
fact, using Dgg = v as predicted by the theory, the match with simulation results over 
the entire range of m is quite good — see Figure [TOl For systems with interacting 
particles, reasonable fits result only if the effective diffusion constant, -Defr, is drastically 
changed from the naively expected value. Moreover, the quality of the the fit is not 
uniformly good over all m. As an example, in Figure [Tl|, where simulation data (for 
L = 10^, i = 10^, u{n) = 1 + 4/n, and p = 0.25, below the condensation threshold) are 
compared with theory plots with three values of D^g, we see that a value of 120 matches 
the data reasonably well. To put this result in context, the naive ("bare") value of the 
diffusion constant is 3.33, so that we would need to invoke "renormalisation effects" at 
the level of a factor of ~ 35. Such a sizable "renormalisation" is comparable to that 
observed in the totally asymmetric simple exclusion process (TASEP) [7j. We believe 
the origin is universal - inter-particle interactions - and that the resolution of this issue 
can be applied to both models. Even treating D^s as a phenomenological parameter, 
the theory fits well only in the midrange, with some discrepancies in the low and high m 
regimes. In the low-m regime the peaks and troughs are not ideally fitted by the theory; 
changing Des to match this region typically destroys the agreement in the mid-range. 
For m > 1500 (not shown here), the lack of agreement is also qualitatively similar, with 
the addition of occasional extra peaks due to a cancellation in the first term in the 
denominator of [T9l 

Since Des is seriously modified by interactions, we investigated systematically the 
effects of varying various parameters in the ZRP. In contrast, changing the parameters 
of the non-interacting system leaves the D^s = v relation completely intact. 

Parameter 6 — The parameter h in the hop rate u = 1 + b/n can be thought of 
as a measure of the strength of the interaction: increasing b first makes a condensation 
transition possible and then reduces pc, the density at which the transition occurs. With 
the density, system size, segment size and normalisation (7) kept constant, the effects 
of varying b are investigated. Now, neither v nor Dcs varies linearly with b, so that 
we find it more meaningful to plot Dcs against v (Figure fT2l) . especially to highlight 
the contrast with D^s = v for a non-interacting system. From the figure, seems to 
change with v mostly in a linear way, although the intercept of this linear component 
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Figure 10. Comparison of power spectra talcen from simulation data and generated 
by the simple theory for the noninteracting system. The system is a segment of £ = 500 
sites in a lattice oi L = 5000 with a global particle density of p — 0.1. Note that in 
the noninteracting system the velocity is 1, but the timescale for the simulation (and 
theory) has been changed such that its actual value is 0.2. For this system the fit is 
apparently good over the entire range and D^s = v. 



— Simulation 




400 800 1200 

m 



Figure 11. Comparison of power spectra taken from simulation data and generated 
by the simple theory for the interacting system with the standard hop rate. The system 
is a segment oi I — 1000 sites in a lattice of L = 10000 with a global particle density 
of 0.25 and hop rate u{n) = 1 + 4/n. The system has a velocity of 3.32653. Here it 
is apparent that the effective diffusion constant is a long way removed from the naive 
expected value but also that the fit is only really good in the midrange of plot, around 
where the second peak due to the size of the segment lies. At the low end, the diffusion 
seems to be too small to capture the rapid oscillations perfectly and at the high end, 
although not shown here, there is an effect from a cancellation in the first term of the 
denominator in p^. 
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Figure 12. Effective diffusion constant plotted against velocity for a varying value of 
the standard hop rate parameter, b. Data shown are fits between the simple theory 
and simulations of a system of i = 10000 sites with a segment of ^ = 1000 at a density 
of p = 0.25 and h varying from 0.5 to 4.0 in increments of 0.5, with low h corresponding 
to low velocity. 



will not go through the origin. Of course, we should expect more serious non-linear 
dependence outside the regime shown here. 

Normalisation 7 — When simulating systems such as the one under consideration 
here, it is common to assign unit probability to the most likely event, in order to 
reduce simulation time. This corresponds to choosing 7 = 1/ maxjw (n)}. For the data 
presented in Figure [12] above, this convention was not followed. Instead, a single 7 was 
used for all the systems (to ensure a conformity of time, t). To explore the effect of 
varying 7, we also carried out simulations with the usual convention. The behaviour 
of Deff with the velocity was found to change dramatically: the previously positive 
trend of Dgfr with v was reversed. To investigate this further, the effect of changing 
7 without other changes was studied. As shown in Figure [131 this gives an apparently 
linear relationship between D^s and v. This is reassuring, since changing 7 alone should 
correspond to nothing more than changing the time scale, while both Defr and v are 
linear in t. More puzzling is the value of the gradient in this plot: The line D^s = 20v 
fits well inside the error bars. Note that a factor of 20 is comparable to the factor 35 
shown earlier. We believe that, once the origin of the substantial renormalisation of Defr 
is uncovered, the resolution of this puzzle will follow. 

Density — With fixed L, i, and b, increasing the overall density (p) lowers both 
V and Deff dramatically. An example of Des vs. p is shown in Figure [Ml (a). Once we 
take into account the v-p relationship and plot Defr against v (Figure [HI b), we again 
recover the line Des ~ 20v, with possibly a small negative curvature. 

Segment Size — Changing the segment size {£) has no effect on the velocity, 
but it does have a pronounced effect on the value of Des, as shown in Figure [151 The 
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Figure 13. Effective diffusion constant plotted against velocity for a varying value of 
the normalisation, 7. Data shown are fits between the simple theory and simulations 
of a system with L = 10000 sites, a segment £ = 1000 sites, at a density of p = 0.1, hop 
rate 1 + 4/n, and 7 values 100-\ 50-\ 25-\ IS'S 13'^, ll-\ 9"^ and 7"^; smaller 
values of 7 correspond to smaller velocities. 



relationship is sub-linear, a behaviour that could be traced to the conserved dynamics. 
That is, the particle numbers in a segment {£) plus those in the complementary segment 
(L — i) is a constant, so that the two averaged power spectra must be the same. 
Unless Dgff (£) develops a singularity at the symmetry point, £ — L/2, we must have 
d^Des {L/2) — 0. Prom this perspective, a sub- linear variation may have been expected 
so as to give a flat profile around i — L/2. However, it is not obvious why the relationship 
should take this form in general. 

System Size — It is also perhaps surprising that varying the system size does 
not appear to change -Dcfr- It is especially so when taken in conjunction with the fact 
that changing the segment size does have an effect. 

It is best to summarise our findings as follows: Although the observed power spectra 
can be reasonably well fitted by the predictions of a linear theory, we must regard 
the effective diffusion constant Des and the noise amplitude A as phenomenological 
parameters. By contrast, the data is entirely consistent with v, the velocity predicted 
from the theory of ZRPs. It is clear that the simulation results do not support the 
prediction of the "naive" theory: Des — v. In addition, by altering the control 
parameters in our study (system size, hopping rate, overall density, and segment size), 
Dcs is not only affected dramatically, but also in such a way that its relationship with 
V changes. At this stage, none of these features are well understood. 
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Figure 14. (a) Effective diffusion constant plotted against density for a varying value 
of the density, p. A non-linear decrease of the effective diffusion constant is shown 
with increasing density, (b) Effective diffusion constant plotted against velocity for a 
varying value of the density, p. Densities are p — 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 
0.6, 0.8, 1.0, 1.6, 2.0, with increasing density corresponding to decreasing velocity. For 
both, data shown are fits between the simple theory and simulations of a system with 
L = 10000 sites, a segment of ^ = 1000 and a hop rate of 1 + 1/n. 



6. Conclusion 



In this paper, we studied the dynamics of fluctuations in the non-equihbrium steady state 
of the zero-range process (ZRP). Speciflcally, we collected time series of the number of 
particles in a contiguous segment of a ring lattice, and computed their average power 
spectra, liuj). We found interesting structures in /(cj), namely, two distinct damped- 
oscillation components. The small uj component consists of narrow peaks over a smooth 
background. The other component resembles broad dips, similar to those observed in 
the power spectra of open TASEP [7]. The origins of these can be traced to, respectively, 



Power Spectra in a ZRP 



19 



300 r 



.1250 
o 

j3200^ 
g 

'gl50^ 



o o 



T " 
o 



I 50^ 



1000 2000 3000 4000 
Segment Size 



5000 



6000 



Figure 15. Effective diffusion constant plotted against segment size. Data stiown are 
fits between tlie simple theory and simulations of a system with L = 10000 sites, a 
density of p = 0.25 and a hop rate of 1 + 4/n. 



the time it takes a fluctuation to travel around the ring and the time for traversing the 
segment. 

We presented a simple toy model, to shed some light on these two types of 
oscillations: a single particle moving ballistically around a ring in continuous space- 
time. The time scries of the "total particle occupation" in a segment is just a periodic 
square wave, so that its Fourier transform is the product of a comb and a sine function, 
controlled respectively by the ring and segment lengths. Diffusion and noise would 
broaden the peaks of the comb and fill in the zeros of the sine function, so that the 
oscillations will appear damped. Let us emphasise that these oscillations are controlled 
by the system and segment sizes, so that they are absent from the usual autocorrelation 
function (for particles at one site, in the thermodynamic limit). Based on the insight 
from the toy model, we believe these features are universal for the power spectra of all 
finite driven diffusive systems. 

At the quantitative level, the observed / (a;) can be fitted quite well by a somewhat 
more sophisticated approach, based on a Langevin equation for the local particle density. 
Focusing on systems far from criticality, we were motivated to linearise this equation 
about the average density. The solution of such an approximate equation, even if 
we account for discrete space-time, is easy. However, except for the case with non- 
interacting particles, not all the parameters of this simple theory fit the simulation 
results. In particular, by identifying the even/odd parity (i.e., V ^ ±V) terms with 
diffusion and drift, we assign the parameters D^s and v, respectively. Good fits can be 
achieved only when D^s is chosen to be considerably larger than the naively predicted 
value. By contrast, v from the simple theory appears to be adequate. At present, we can 
only present the dependence of D^s on various control parameters as phenomenological 
results from our extensive simulation studies. In the same vein, the relationship between 
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Deflf and V (for the range of systems we considered) was seen to deviate significantly 
from the naive theoretical prediction of Defr = v. Such "anomalies" associated with D^^ 
were also observed in similar studies of another system [7j. Another puzzling aspect 
here is that Des appears to depend more strongly on the segment size {€) than on the 
system size (L). This feature may indicate that the introduction of an effective diffusion 
constant to a linear Langevin theory is not an entirely satisfactory treatment. While 
this approach is somewhat successful at fitting individual power spectra, it leaves much 
room for improvement, as we seek better understanding and a comprehensive theory. 

Although we have ruled out the possibility that these anomalies are due to a 
systematic effect of our random number generator (by using different generators and 
choosing parameters to be incommensurate with each other), we have not considered 
alternative simulation methods. Two such alternatives come readily to mind. One is 
kinetic Monte Carlo [18], where an appropriate event is chosen at each update and 
the time advanced according to a Poisson distribution. Another alternative would 
be to try a common method of reducing simulation times which is to pick a particle 
at random and hop with an appropriate probability. Fundamentally, we believe that 
inter-particle interactions are responsible for large deviations from the simple linearised 
theory presented here, rather than some subtle effect due to the details of the particular 
dynamics we chose. 

It is interesting that similar structures in the power spectra have been observed 
in two of the most simple models for non-equilibrium systems, namely TASEP and 
ZRP. Oscillations seen in the variance of the integrated current at a single site in the 
time domain [10] are also undoubtedly related. These suggest that damped-oscillatory 
behaviour is universal for finite systems driven out of equilibrium. Further investigations 
to place this notion on a sound foundation would be worthwhile. It is clear that 
fluctuations in non-equilibrium steady states are non-trivial and their dynamics induce 
interesting behaviour. This study has shown only a limited view in a small corner of 
this vast area. Even within this corner, there is room for improvements, especially in 
more complete analytic theory. We hope that a better understanding of this particular 
problem will lead to deeper insights into the nature of fluctuations in physical systems 
driven far from thermal equilibrium. 
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